Polynommultiplikation

polymult <- function(x,a,b)
  # a(B) = a_0 + a_1*B + a_2*B² + ... + a_p*B^p
  # b(B) = b_0 + b_1*B + b_2*B² + ... + a_q*B^q
  # where wlog p >= q
  # Result is a(B)b(B) x_t
{
  if (length(b) == 1){polycoeff <- b[1]*array(a)}
  else
  {alpha <- array(a)
  p <- dim(alpha) - 1
  beta <- array(b)
  q <- dim(beta) - 1
  # Assumed: p >= q
  polycoeff <- array(rep(0,p+q + 1))
  for (k in 0:(p+q))
    for (j in max(0,k-q):min(k,p))
    {polycoeff[k+1] <- polycoeff[k+1] + alpha[j+1]*beta[k-j+1]}
  }
  return(filter(x,c(polycoeff),sides=1, method="convolution"))
  # test: return(polycoeff)
}

Opstilling af Data

library(readxl)
library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
library(imputeTS)
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'imputeTS'
## The following object is masked from 'package:zoo':
## 
##     na.locf
library(forecast)
## Warning: package 'forecast' was built under R version 3.6.2
library(tseries)
## 
## Attaching package: 'tseries'
## The following object is masked from 'package:imputeTS':
## 
##     na.remove
library(stats)
library(strucchange)
library(tidyr)
library(astsa)
## 
## Attaching package: 'astsa'
## The following object is masked from 'package:forecast':
## 
##     gas
elprice17 =  data.frame(read.csv("elspot-prices_2017_daily_eur.csv",sep=";",dec=","))

cons17= data.frame(read.csv("consumption-se-areas_2017_daily.csv",sep=";"))

prod17 = data.frame(read.csv("production-se-areas_2017_daily.csv",sep=";"))


Price = log(na_interpolation(as.numeric(elprice17[3:367,4]), option = "linear"))

Cons = log(na_interpolation(as.numeric(as.character(cons17[3:367,3])), option = "linear"))

Prod = log(na_interpolation(as.numeric(as.character(prod17[3:367,3])), option = "linear"))


plot.ts(Price)

plot.ts(Cons)

plot.ts(Prod)

acf(Price)

acf(Cons)

acf(Prod)

dato <- seq(c(ISOdate(2017,1,1)), by = "days", length.out = 365)

data1 = data.frame(dato, Price, Cons, Prod)

Prod1= data1[,c(1,4)]

Cons1 = data1[,c(1,3)]

Price1 = data1[,c(1,2)]

Estimation af deterministisk sæson

sprod = glm(Prod1[,2] ~ time(Prod1[,1]) + 
                    I(time(Prod1[,1])^2) +
                    sin((4*pi)/365*I(time(Prod1[,1])))+ 
                    cos((4*pi)/365*I(time(Prod1[,1])))+  
                    sin((8*pi)/365*I(time(Prod1[,1])))+ 
                    cos((8*pi)/365*I(time(Prod1[,1])))+
                    sin((24*pi)/365*I(time(Prod1[,1])))+ 
                    cos((24*pi)/365*I(time(Prod1[,1])))+
                    sin((104*pi)/365*I(time(Prod1[,1])))+ 
                    cos((104*pi)/365*I(time(Prod1[,1])))+
                    sin((730*pi)/365*I(time(Prod1[,1])))+ 
                    cos((730*pi)/365*I(time(Prod1[,1])))
)

sprice = glm(Price1[,2] ~ time(Price1[,1]) + 
              I(time(Price1[,1])^2) +
              sin((4*pi)/365*I(time(Price1[,1])))+ 
              cos((4*pi)/365*I(time(Price1[,1])))+  
              sin((8*pi)/365*I(time(Price1[,1])))+ 
              cos((8*pi)/365*I(time(Price1[,1])))+
              sin((24*pi)/365*I(time(Price1[,1])))+ 
              cos((24*pi)/365*I(time(Price1[,1])))+
              sin((104*pi)/365*I(time(Price1[,1])))+ 
              cos((104*pi)/365*I(time(Price1[,1])))+
              sin((730*pi)/365*I(time(Price1[,1])))+ 
              cos((730*pi)/365*I(time(Price1[,1])))
)


scons = glm(Cons1[,2] ~ time(Cons1[,1]) + 
              I(time(Cons1[,1])^2) +
              sin((4*pi)/365*I(time(Cons1[,1])))+ 
              cos((4*pi)/365*I(time(Cons1[,1])))+  
              sin((8*pi)/365*I(time(Cons1[,1])))+ 
              cos((8*pi)/365*I(time(Cons1[,1])))+
              sin((24*pi)/365*I(time(Cons1[,1])))+ 
              cos((24*pi)/365*I(time(Cons1[,1])))+
              sin((104*pi)/365*I(time(Cons1[,1])))+ 
              cos((104*pi)/365*I(time(Cons1[,1])))+
              sin((730*pi)/365*I(time(Cons1[,1])))+ 
              cos((730*pi)/365*I(time(Cons1[,1])))
)

summary(sprod)
## 
## Call:
## glm(formula = Prod1[, 2] ~ time(Prod1[, 1]) + I(time(Prod1[, 
##     1])^2) + sin((4 * pi)/365 * I(time(Prod1[, 1]))) + cos((4 * 
##     pi)/365 * I(time(Prod1[, 1]))) + sin((8 * pi)/365 * I(time(Prod1[, 
##     1]))) + cos((8 * pi)/365 * I(time(Prod1[, 1]))) + sin((24 * 
##     pi)/365 * I(time(Prod1[, 1]))) + cos((24 * pi)/365 * I(time(Prod1[, 
##     1]))) + sin((104 * pi)/365 * I(time(Prod1[, 1]))) + cos((104 * 
##     pi)/365 * I(time(Prod1[, 1]))) + sin((730 * pi)/365 * I(time(Prod1[, 
##     1]))) + cos((730 * pi)/365 * I(time(Prod1[, 1]))))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.44968  -0.08051   0.01333   0.08967   0.32653  
## 
## Coefficients: (1 not defined because of singularities)
##                                             Estimate Std. Error t value
## (Intercept)                                1.187e+01  2.253e-02 527.138
## time(Prod1[, 1])                          -3.242e-03  2.879e-04 -11.262
## I(time(Prod1[, 1])^2)                      8.745e-06  7.346e-07  11.905
## sin((4 * pi)/365 * I(time(Prod1[, 1])))    4.048e-02  1.090e-02   3.715
## cos((4 * pi)/365 * I(time(Prod1[, 1])))   -1.521e-02  1.030e-02  -1.477
## sin((8 * pi)/365 * I(time(Prod1[, 1])))    1.824e-02  1.023e-02   1.784
## cos((8 * pi)/365 * I(time(Prod1[, 1])))   -2.846e-02  1.001e-02  -2.842
## sin((24 * pi)/365 * I(time(Prod1[, 1])))  -5.563e-03  1.002e-02  -0.555
## cos((24 * pi)/365 * I(time(Prod1[, 1])))  -3.212e-02  9.995e-03  -3.213
## sin((104 * pi)/365 * I(time(Prod1[, 1]))) -8.109e-03  1.000e-02  -0.811
## cos((104 * pi)/365 * I(time(Prod1[, 1]))) -1.009e-01  9.996e-03 -10.097
## sin((730 * pi)/365 * I(time(Prod1[, 1])))  7.334e+09  1.106e+11   0.066
## cos((730 * pi)/365 * I(time(Prod1[, 1])))         NA         NA      NA
##                                           Pr(>|t|)    
## (Intercept)                                < 2e-16 ***
## time(Prod1[, 1])                           < 2e-16 ***
## I(time(Prod1[, 1])^2)                      < 2e-16 ***
## sin((4 * pi)/365 * I(time(Prod1[, 1])))   0.000236 ***
## cos((4 * pi)/365 * I(time(Prod1[, 1])))   0.140549    
## sin((8 * pi)/365 * I(time(Prod1[, 1])))   0.075320 .  
## cos((8 * pi)/365 * I(time(Prod1[, 1])))   0.004739 ** 
## sin((24 * pi)/365 * I(time(Prod1[, 1])))  0.579141    
## cos((24 * pi)/365 * I(time(Prod1[, 1])))  0.001432 ** 
## sin((104 * pi)/365 * I(time(Prod1[, 1]))) 0.417974    
## cos((104 * pi)/365 * I(time(Prod1[, 1])))  < 2e-16 ***
## sin((730 * pi)/365 * I(time(Prod1[, 1]))) 0.947152    
## cos((730 * pi)/365 * I(time(Prod1[, 1])))       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.01822984)
## 
##     Null deviance: 11.6055  on 364  degrees of freedom
## Residual deviance:  6.4351  on 353  degrees of freedom
## AIC: -412.09
## 
## Number of Fisher Scoring iterations: 2
summary(scons)
## 
## Call:
## glm(formula = Cons1[, 2] ~ time(Cons1[, 1]) + I(time(Cons1[, 
##     1])^2) + sin((4 * pi)/365 * I(time(Cons1[, 1]))) + cos((4 * 
##     pi)/365 * I(time(Cons1[, 1]))) + sin((8 * pi)/365 * I(time(Cons1[, 
##     1]))) + cos((8 * pi)/365 * I(time(Cons1[, 1]))) + sin((24 * 
##     pi)/365 * I(time(Cons1[, 1]))) + cos((24 * pi)/365 * I(time(Cons1[, 
##     1]))) + sin((104 * pi)/365 * I(time(Cons1[, 1]))) + cos((104 * 
##     pi)/365 * I(time(Cons1[, 1]))) + sin((730 * pi)/365 * I(time(Cons1[, 
##     1]))) + cos((730 * pi)/365 * I(time(Cons1[, 1]))))
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.204789  -0.043860   0.003861   0.044519   0.236736  
## 
## Coefficients: (1 not defined because of singularities)
##                                             Estimate Std. Error  t value
## (Intercept)                                1.123e+01  1.119e-02 1004.281
## time(Cons1[, 1])                          -7.021e-03  1.430e-04  -49.108
## I(time(Cons1[, 1])^2)                      1.691e-05  3.648e-07   46.347
## sin((4 * pi)/365 * I(time(Cons1[, 1])))   -2.616e-02  5.411e-03   -4.834
## cos((4 * pi)/365 * I(time(Cons1[, 1])))   -2.665e-02  5.114e-03   -5.211
## sin((8 * pi)/365 * I(time(Cons1[, 1])))   -1.762e-02  5.079e-03   -3.470
## cos((8 * pi)/365 * I(time(Cons1[, 1])))   -9.944e-03  4.973e-03   -2.000
## sin((24 * pi)/365 * I(time(Cons1[, 1])))  -8.262e-04  4.976e-03   -0.166
## cos((24 * pi)/365 * I(time(Cons1[, 1])))  -6.038e-03  4.964e-03   -1.216
## sin((104 * pi)/365 * I(time(Cons1[, 1]))) -2.717e-03  4.966e-03   -0.547
## cos((104 * pi)/365 * I(time(Cons1[, 1]))) -4.006e-02  4.964e-03   -8.070
## sin((730 * pi)/365 * I(time(Cons1[, 1])))  1.251e+11  5.491e+10    2.278
## cos((730 * pi)/365 * I(time(Cons1[, 1])))         NA         NA       NA
##                                           Pr(>|t|)    
## (Intercept)                                < 2e-16 ***
## time(Cons1[, 1])                           < 2e-16 ***
## I(time(Cons1[, 1])^2)                      < 2e-16 ***
## sin((4 * pi)/365 * I(time(Cons1[, 1])))   2.00e-06 ***
## cos((4 * pi)/365 * I(time(Cons1[, 1])))   3.21e-07 ***
## sin((8 * pi)/365 * I(time(Cons1[, 1])))   0.000585 ***
## cos((8 * pi)/365 * I(time(Cons1[, 1])))   0.046306 *  
## sin((24 * pi)/365 * I(time(Cons1[, 1])))  0.868234    
## cos((24 * pi)/365 * I(time(Cons1[, 1])))  0.224641    
## sin((104 * pi)/365 * I(time(Cons1[, 1]))) 0.584693    
## cos((104 * pi)/365 * I(time(Cons1[, 1]))) 1.11e-14 ***
## sin((730 * pi)/365 * I(time(Cons1[, 1]))) 0.023312 *  
## cos((730 * pi)/365 * I(time(Cons1[, 1])))       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.004495528)
## 
##     Null deviance: 13.6568  on 364  degrees of freedom
## Residual deviance:  1.5869  on 353  degrees of freedom
## AIC: -923.08
## 
## Number of Fisher Scoring iterations: 2
summary(sprice)
## 
## Call:
## glm(formula = Price1[, 2] ~ time(Price1[, 1]) + I(time(Price1[, 
##     1])^2) + sin((4 * pi)/365 * I(time(Price1[, 1]))) + cos((4 * 
##     pi)/365 * I(time(Price1[, 1]))) + sin((8 * pi)/365 * I(time(Price1[, 
##     1]))) + cos((8 * pi)/365 * I(time(Price1[, 1]))) + sin((24 * 
##     pi)/365 * I(time(Price1[, 1]))) + cos((24 * pi)/365 * I(time(Price1[, 
##     1]))) + sin((104 * pi)/365 * I(time(Price1[, 1]))) + cos((104 * 
##     pi)/365 * I(time(Price1[, 1]))) + sin((730 * pi)/365 * I(time(Price1[, 
##     1]))) + cos((730 * pi)/365 * I(time(Price1[, 1]))))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4750  -0.4025   0.1499   0.5555   1.3704  
## 
## Coefficients: (1 not defined because of singularities)
##                                              Estimate Std. Error t value
## (Intercept)                                 5.063e+00  1.371e-01  36.936
## time(Price1[, 1])                          -5.554e-03  1.752e-03  -3.170
## I(time(Price1[, 1])^2)                      1.726e-05  4.471e-06   3.860
## sin((4 * pi)/365 * I(time(Price1[, 1])))    3.042e-01  6.631e-02   4.588
## cos((4 * pi)/365 * I(time(Price1[, 1])))   -1.958e-01  6.267e-02  -3.125
## sin((8 * pi)/365 * I(time(Price1[, 1])))   -1.405e-02  6.224e-02  -0.226
## cos((8 * pi)/365 * I(time(Price1[, 1])))   -1.375e-01  6.094e-02  -2.256
## sin((24 * pi)/365 * I(time(Price1[, 1])))   8.059e-02  6.098e-02   1.322
## cos((24 * pi)/365 * I(time(Price1[, 1])))  -1.240e-02  6.083e-02  -0.204
## sin((104 * pi)/365 * I(time(Price1[, 1]))) -3.065e-02  6.086e-02  -0.504
## cos((104 * pi)/365 * I(time(Price1[, 1]))) -3.848e-01  6.083e-02  -6.326
## sin((730 * pi)/365 * I(time(Price1[, 1]))) -9.109e+10  6.729e+11  -0.135
## cos((730 * pi)/365 * I(time(Price1[, 1])))         NA         NA      NA
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## time(Price1[, 1])                          0.001656 ** 
## I(time(Price1[, 1])^2)                     0.000135 ***
## sin((4 * pi)/365 * I(time(Price1[, 1])))   6.22e-06 ***
## cos((4 * pi)/365 * I(time(Price1[, 1])))   0.001928 ** 
## sin((8 * pi)/365 * I(time(Price1[, 1])))   0.821594    
## cos((8 * pi)/365 * I(time(Price1[, 1])))   0.024676 *  
## sin((24 * pi)/365 * I(time(Price1[, 1])))  0.187161    
## cos((24 * pi)/365 * I(time(Price1[, 1])))  0.838569    
## sin((104 * pi)/365 * I(time(Price1[, 1]))) 0.614818    
## cos((104 * pi)/365 * I(time(Price1[, 1]))) 7.64e-10 ***
## sin((730 * pi)/365 * I(time(Price1[, 1]))) 0.892398    
## cos((730 * pi)/365 * I(time(Price1[, 1])))       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.6751302)
## 
##     Null deviance: 297.25  on 364  degrees of freedom
## Residual deviance: 238.32  on 353  degrees of freedom
## AIC: 906.23
## 
## Number of Fisher Scoring iterations: 2
sprod_res = as.ts(sprod$residuals)
sprice_res = as.ts(sprice$residuals)
scons_res = as.ts(scons$residuals)

acf2(sprod_res)

##         ACF  PACF
##  [1,]  0.49  0.49
##  [2,]  0.24  0.01
##  [3,]  0.31  0.24
##  [4,]  0.23 -0.01
##  [5,]  0.06 -0.10
##  [6,]  0.12  0.11
##  [7,]  0.30  0.24
##  [8,]  0.05 -0.26
##  [9,] -0.03  0.02
## [10,]  0.12  0.07
## [11,]  0.10  0.02
## [12,] -0.03 -0.04
## [13,]  0.03  0.02
## [14,]  0.21  0.13
## [15,] -0.03 -0.18
## [16,] -0.09 -0.01
## [17,]  0.05  0.02
## [18,] -0.02 -0.08
## [19,] -0.14 -0.03
## [20,] -0.13 -0.12
## [21,] -0.01  0.02
## [22,] -0.23 -0.17
## [23,] -0.31 -0.10
## [24,] -0.19 -0.07
## [25,] -0.17  0.04
## [26,] -0.27 -0.08
## [27,] -0.21 -0.02
## [28,] -0.06  0.00
## [29,] -0.26 -0.14
## [30,] -0.33 -0.08
acf2(sprice_res)

##         ACF  PACF
##  [1,]  0.53  0.53
##  [2,]  0.17 -0.16
##  [3,]  0.16  0.20
##  [4,]  0.12 -0.06
##  [5,]  0.04  0.00
##  [6,]  0.11  0.14
##  [7,]  0.20  0.09
##  [8,]  0.04 -0.16
##  [9,] -0.01  0.08
## [10,]  0.09  0.06
## [11,]  0.08 -0.02
## [12,]  0.01 -0.01
## [13,]  0.04  0.03
## [14,]  0.07 -0.01
## [15,] -0.04 -0.06
## [16,] -0.02  0.05
## [17,]  0.03 -0.04
## [18,] -0.04 -0.06
## [19,] -0.10 -0.03
## [20,] -0.04  0.02
## [21,] -0.02 -0.04
## [22,] -0.14 -0.10
## [23,] -0.15 -0.04
## [24,] -0.13 -0.08
## [25,] -0.13  0.03
## [26,] -0.09  0.01
## [27,] -0.07 -0.07
## [28,] -0.10 -0.03
## [29,] -0.15 -0.02
## [30,] -0.16 -0.09
acf2(scons_res)

##         ACF  PACF
##  [1,]  0.61  0.61
##  [2,]  0.40  0.04
##  [3,]  0.31  0.07
##  [4,]  0.20 -0.03
##  [5,]  0.13 -0.01
##  [6,]  0.17  0.13
##  [7,]  0.26  0.18
##  [8,]  0.21 -0.07
##  [9,]  0.19  0.03
## [10,]  0.19  0.03
## [11,]  0.12 -0.04
## [12,]  0.00 -0.12
## [13,]  0.00  0.03
## [14,]  0.05  0.06
## [15,]  0.00 -0.07
## [16,]  0.02  0.02
## [17,]  0.03 -0.05
## [18,]  0.01  0.01
## [19,] -0.02  0.01
## [20,] -0.03 -0.03
## [21,]  0.01  0.05
## [22,] -0.01 -0.01
## [23,]  0.01  0.04
## [24,]  0.09  0.09
## [25,]  0.08 -0.01
## [26,]  0.04 -0.03
## [27,]  0.00 -0.06
## [28,]  0.03  0.08
## [29,] -0.06 -0.15
## [30,] -0.08 -0.02

Estimation og valg af SARMA for Pris

acf2(sprice_res)

##         ACF  PACF
##  [1,]  0.53  0.53
##  [2,]  0.17 -0.16
##  [3,]  0.16  0.20
##  [4,]  0.12 -0.06
##  [5,]  0.04  0.00
##  [6,]  0.11  0.14
##  [7,]  0.20  0.09
##  [8,]  0.04 -0.16
##  [9,] -0.01  0.08
## [10,]  0.09  0.06
## [11,]  0.08 -0.02
## [12,]  0.01 -0.01
## [13,]  0.04  0.03
## [14,]  0.07 -0.01
## [15,] -0.04 -0.06
## [16,] -0.02  0.05
## [17,]  0.03 -0.04
## [18,] -0.04 -0.06
## [19,] -0.10 -0.03
## [20,] -0.04  0.02
## [21,] -0.02 -0.04
## [22,] -0.14 -0.10
## [23,] -0.15 -0.04
## [24,] -0.13 -0.08
## [25,] -0.13  0.03
## [26,] -0.09  0.01
## [27,] -0.07 -0.07
## [28,] -0.10 -0.03
## [29,] -0.15 -0.02
## [30,] -0.16 -0.09
seas_sprice_res1 = sarima(sprice_res,3,0,0,1,0,0,7, no.constant = TRUE);
## initial  value -0.205739 
## iter   2 value -0.329456
## iter   3 value -0.412408
## iter   4 value -0.422958
## iter   5 value -0.425715
## iter   6 value -0.425722
## iter   7 value -0.425722
## iter   7 value -0.425722
## iter   7 value -0.425722
## final  value -0.425722 
## converged
## initial  value -0.433484 
## iter   2 value -0.433492
## iter   3 value -0.433492
## iter   3 value -0.433492
## iter   3 value -0.433492
## final  value -0.433492 
## converged

seas_sprice_res2 = sarima(sprice_res,2,0,0,1,0,0,7, no.constant = TRUE);
## initial  value -0.206984 
## iter   2 value -0.355555
## iter   3 value -0.409459
## iter   4 value -0.414653
## iter   5 value -0.414788
## iter   6 value -0.414788
## iter   7 value -0.414788
## iter   7 value -0.414788
## iter   7 value -0.414788
## final  value -0.414788 
## converged
## initial  value -0.421343 
## iter   2 value -0.421358
## iter   3 value -0.421358
## iter   3 value -0.421358
## iter   3 value -0.421358
## final  value -0.421358 
## converged

seas_sprice_res3 = sarima(sprice_res,4,0,0,1,0,0,7, no.constant = TRUE);
## initial  value -0.206440 
## iter   2 value -0.305208
## iter   3 value -0.409103
## iter   4 value -0.421271
## iter   5 value -0.427345
## iter   6 value -0.428170
## iter   7 value -0.428190
## iter   8 value -0.428190
## iter   8 value -0.428190
## iter   8 value -0.428190
## final  value -0.428190 
## converged
## initial  value -0.435278 
## iter   2 value -0.435301
## iter   3 value -0.435301
## iter   4 value -0.435301
## iter   4 value -0.435301
## iter   4 value -0.435301
## final  value -0.435301 
## converged

seas_sprice_res4 = sarima(sprice_res,3,0,0,2,0,0,7, no.constant = TRUE);
## initial  value -0.203956 
## iter   2 value -0.326430
## iter   3 value -0.412469
## iter   4 value -0.423640
## iter   5 value -0.426529
## iter   6 value -0.426565
## iter   7 value -0.426567
## iter   8 value -0.426567
## iter   8 value -0.426567
## iter   8 value -0.426567
## final  value -0.426567 
## converged
## initial  value -0.435801 
## iter   2 value -0.435814
## iter   3 value -0.435814
## iter   4 value -0.435814
## iter   4 value -0.435814
## iter   4 value -0.435814
## final  value -0.435814 
## converged

seas_sprice_res5 = sarima(sprice_res,3,0,0,0,0,0,7, no.constant = TRUE);
## initial  value -0.210262 
## iter   2 value -0.312482
## iter   3 value -0.391609
## iter   4 value -0.406578
## iter   5 value -0.408806
## iter   6 value -0.408817
## iter   6 value -0.408817
## iter   6 value -0.408817
## final  value -0.408817 
## converged
## initial  value -0.410536 
## iter   2 value -0.410539
## iter   3 value -0.410540
## iter   4 value -0.410540
## iter   4 value -0.410540
## iter   4 value -0.410540
## final  value -0.410540 
## converged

seas_sprice_res1$AIC
## [1] 0.1526423
seas_sprice_res2$AIC
## [1] 0.1713686
seas_sprice_res3$AIC
## [1] 0.1544827
seas_sprice_res4$AIC
## [1] 0.1532696
seas_sprice_res5$AIC
## [1] 0.1939803

Estimation og valg af SARMA for Forbrug

seas_scons_res1 = sarima(scons_res,1,0,0,1,0,0,7, no.constant = TRUE)
## initial  value -2.731571 
## iter   2 value -3.000725
## iter   3 value -3.002720
## iter   4 value -3.002761
## iter   5 value -3.002772
## iter   5 value -3.002772
## iter   5 value -3.002772
## final  value -3.002772 
## converged
## initial  value -2.987264 
## iter   2 value -2.987335
## iter   3 value -2.987335
## iter   3 value -2.987335
## iter   3 value -2.987335
## final  value -2.987335 
## converged

seas_scons_res2 = sarima(scons_res,2,0,0,1,0,0,7, no.constant = TRUE)
## initial  value -2.732810 
## iter   2 value -2.878819
## iter   3 value -2.992100
## iter   4 value -3.001923
## iter   5 value -3.003134
## iter   6 value -3.003137
## iter   7 value -3.003137
## iter   7 value -3.003137
## iter   7 value -3.003137
## final  value -3.003137 
## converged
## initial  value -2.987268 
## iter   2 value -2.987405
## iter   3 value -2.987408
## iter   4 value -2.987408
## iter   4 value -2.987408
## iter   4 value -2.987408
## final  value -2.987408 
## converged

seas_scons_res3 = sarima(scons_res,1,0,0,2,0,0,7, no.constant = TRUE)
## initial  value -2.737468 
## iter   2 value -3.005414
## iter   3 value -3.008169
## iter   4 value -3.008307
## iter   5 value -3.008308
## iter   6 value -3.008309
## iter   6 value -3.008309
## iter   6 value -3.008309
## final  value -3.008309 
## converged
## initial  value -2.990452 
## iter   2 value -2.990510
## iter   3 value -2.990513
## iter   4 value -2.990514
## iter   4 value -2.990514
## iter   4 value -2.990514
## final  value -2.990514 
## converged

seas_scons_res4 = sarima(scons_res,1,0,0,0,0,0,7, no.constant = TRUE)
## initial  value -2.731071 
## iter   2 value -2.976754
## iter   3 value -2.976883
## iter   4 value -2.976959
## iter   4 value -2.976959
## final  value -2.976959 
## converged
## initial  value -2.964079 
## iter   2 value -2.964230
## iter   3 value -2.964257
## iter   4 value -2.964257
## iter   4 value -2.964257
## final  value -2.964257 
## converged

seas_scons_res1$AIC
## [1] -4.965997
seas_scons_res2$AIC 
## [1] -4.960661
seas_scons_res3$AIC 
## [1] -4.96718
seas_scons_res4$AIC 
## [1] -4.924423

Estimation og valg af SARMA for Produktion

seas_sprod_res1 = sarima(sprod_res,1,0,0,2,0,0,7, no.constant = TRUE)
## initial  value -2.031458 
## iter   2 value -2.303345
## iter   3 value -2.304909
## iter   4 value -2.310272
## iter   5 value -2.310292
## iter   6 value -2.310301
## iter   7 value -2.310301
## iter   7 value -2.310301
## iter   7 value -2.310301
## final  value -2.310301 
## converged
## initial  value -2.280917 
## iter   2 value -2.280944
## iter   3 value -2.280945
## iter   4 value -2.280945
## iter   4 value -2.280945
## iter   4 value -2.280945
## final  value -2.280945 
## converged

seas_sprod_res2 = sarima(sprod_res,2,0,0,2,0,0,7, no.constant = TRUE)
## initial  value -2.031145 
## iter   2 value -2.270979
## iter   3 value -2.297740
## iter   4 value -2.308261
## iter   5 value -2.308848
## iter   6 value -2.309075
## iter   7 value -2.309075
## iter   7 value -2.309075
## iter   7 value -2.309075
## final  value -2.309075 
## converged
## initial  value -2.281094 
## iter   2 value -2.281176
## iter   3 value -2.281180
## iter   4 value -2.281180
## iter   4 value -2.281180
## iter   4 value -2.281180
## final  value -2.281180 
## converged

seas_sprod_res3 = sarima(sprod_res,3,0,0,2,0,0,7, no.constant = TRUE)
## initial  value -2.029721 
## iter   2 value -2.184333
## iter   3 value -2.303575
## iter   4 value -2.315760
## iter   5 value -2.317568
## iter   6 value -2.318229
## iter   7 value -2.318261
## iter   8 value -2.318261
## iter   8 value -2.318261
## iter   8 value -2.318261
## final  value -2.318261 
## converged
## initial  value -2.289435 
## iter   2 value -2.289565
## iter   3 value -2.289573
## iter   4 value -2.289575
## iter   5 value -2.289575
## iter   5 value -2.289575
## iter   5 value -2.289575
## final  value -2.289575 
## converged

seas_sprod_res4 = sarima(sprod_res,0,0,0,2,0,0,7, no.constant = TRUE)
## initial  value -2.031748 
## iter   2 value -2.088725
## iter   3 value -2.094309
## iter   4 value -2.094449
## iter   5 value -2.094457
## iter   5 value -2.094457
## iter   5 value -2.094457
## final  value -2.094457 
## converged
## initial  value -2.079786 
## iter   2 value -2.079795
## iter   3 value -2.079795
## iter   3 value -2.079795
## iter   3 value -2.079795
## final  value -2.079795 
## converged

seas_sprod_res5 = sarima(sprod_res,1,0,0,1,0,0,7, no.constant = TRUE)
## initial  value -2.032823 
## iter   2 value -2.265262
## iter   3 value -2.271592
## iter   4 value -2.273270
## iter   5 value -2.273276
## iter   6 value -2.273277
## iter   6 value -2.273277
## final  value -2.273277 
## converged
## initial  value -2.250554 
## iter   2 value -2.250598
## iter   3 value -2.250600
## iter   4 value -2.250600
## iter   4 value -2.250600
## iter   4 value -2.250600
## final  value -2.250600 
## converged

seas_sprod_res6 = sarima(sprod_res,1,0,0,0,0,0,7, no.constant = TRUE)
## initial  value -2.033654 
## iter   2 value -2.174635
## iter   3 value -2.174745
## iter   4 value -2.174757
## iter   4 value -2.174757
## final  value -2.174757 
## converged
## initial  value -2.159654 
## iter   2 value -2.159780
## iter   3 value -2.159789
## iter   3 value -2.159789
## iter   3 value -2.159789
## final  value -2.159789 
## converged

seas_sprod_res7 = sarima(sprod_res,1,0,0,3,0,0,7, no.constant = TRUE)
## initial  value -2.026230 
## iter   2 value -2.297344
## iter   3 value -2.301438
## iter   4 value -2.306807
## iter   5 value -2.306896
## iter   6 value -2.306903
## iter   7 value -2.306918
## iter   7 value -2.306918
## iter   7 value -2.306918
## final  value -2.306918 
## converged
## initial  value -2.282058 
## iter   2 value -2.282158
## iter   3 value -2.282176
## iter   4 value -2.282180
## iter   5 value -2.282180
## iter   5 value -2.282180
## iter   5 value -2.282180
## final  value -2.282180 
## converged

seas_sprod_res1$AIC
## [1] -3.553287
seas_sprod_res2$AIC
## [1] -3.548292
seas_sprod_res3$AIC
## [1] -3.559168
seas_sprod_res4$AIC
## [1] -3.151538
seas_sprod_res5$AIC
## [1] -3.49516
seas_sprod_res6$AIC
## [1] -3.314901
seas_sprod_res7$AIC
## [1] -3.550507

SARMA

seas_y1res = ts(seas_sprice_res1$fit$residuals)
coefs_y1 = seas_sprice_res1$fit$coef

coefs_season_Phi_y1 = c(1,rep(0,6),-coefs_y1[2])

stat_y1 = polymult(sprice_res,coefs_season_Phi_y1,1)


seas_x1res = ts(seas_sprod_res3$fit$residuals)
coefs_x1 = seas_sprod_res3$fit$coef

coefs_season_Phi_x1 = c(1,rep(0,6),-coefs_x1[2])

stat_x1 = polymult(sprod_res,coefs_season_Phi_x1,1)


seas_z1res = ts(seas_scons_res3$fit$residuals)
coefs_z1 = seas_scons_res3$fit$coef

coefs_season_Phi_z1 = c(1,rep(0,6),-coefs_z1[2])

stat_z1 = polymult(scons_res,coefs_season_Phi_z1,1)

Opstilling af VAR-model

XXX = data.frame(stat_z1,stat_x1,stat_y1)[8:365,]

var_xx = VAR(XXX, lag.max = 8)

VARselect(XXX)
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      8      8      1      8 
## 
## $criteria
##                    1             2             3             4
## AIC(n) -1.130626e+01 -1.131529e+01 -1.140246e+01 -1.138376e+01
## HQ(n)  -1.125338e+01 -1.122274e+01 -1.127025e+01 -1.121188e+01
## SC(n)  -1.117343e+01 -1.108283e+01 -1.107038e+01 -1.095204e+01
## FPE(n)  1.229569e-05  1.218541e-05  1.116846e-05  1.138002e-05
##                    5             6             7             8
## AIC(n) -1.134625e+01 -1.138470e+01 -1.151625e+01 -1.171132e+01
## HQ(n)  -1.113471e+01 -1.113351e+01 -1.122539e+01 -1.138080e+01
## SC(n)  -1.081491e+01 -1.075374e+01 -1.078566e+01 -1.088111e+01
## FPE(n)  1.181603e-05  1.137175e-05  9.971881e-06  8.206550e-06
##                    9            10
## AIC(n) -1.168681e+01 -1.165476e+01
## HQ(n)  -1.131662e+01 -1.124491e+01
## SC(n)  -1.075697e+01 -1.062529e+01
## FPE(n)  8.412736e-06  8.689992e-06
summary(var_xx)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: stat_z1, stat_x1, stat_y1 
## Deterministic variables: const 
## Sample size: 350 
## Log Likelihood: 636.669 
## Roots of the characteristic polynomial:
## 0.9503 0.9503 0.9012 0.9012 0.8972 0.8972 0.8777 0.8777 0.8729 0.8729 0.8428 0.8428 0.8148 0.8148 0.8133 0.8133 0.8107 0.8107 0.8031 0.782 0.782 0.7655 0.7655 0.7362
## Call:
## VAR(y = XXX, lag.max = 8)
## 
## 
## Estimation results for equation stat_z1: 
## ======================================== 
## stat_z1 = stat_z1.l1 + stat_x1.l1 + stat_y1.l1 + stat_z1.l2 + stat_x1.l2 + stat_y1.l2 + stat_z1.l3 + stat_x1.l3 + stat_y1.l3 + stat_z1.l4 + stat_x1.l4 + stat_y1.l4 + stat_z1.l5 + stat_x1.l5 + stat_y1.l5 + stat_z1.l6 + stat_x1.l6 + stat_y1.l6 + stat_z1.l7 + stat_x1.l7 + stat_y1.l7 + stat_z1.l8 + stat_x1.l8 + stat_y1.l8 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## stat_z1.l1  0.5951356  0.0573602  10.375  < 2e-16 ***
## stat_x1.l1  0.0196595  0.0279328   0.704  0.48205    
## stat_y1.l1 -0.0039777  0.0044897  -0.886  0.37629    
## stat_z1.l2  0.0760000  0.0662359   1.147  0.25205    
## stat_x1.l2 -0.0131725  0.0298994  -0.441  0.65982    
## stat_y1.l2 -0.0066493  0.0051514  -1.291  0.19770    
## stat_z1.l3 -0.0575799  0.0670678  -0.859  0.39123    
## stat_x1.l3  0.0284036  0.0299568   0.948  0.34376    
## stat_y1.l3  0.0075568  0.0052586   1.437  0.15167    
## stat_z1.l4 -0.0084991  0.0670925  -0.127  0.89927    
## stat_x1.l4  0.0234458  0.0300899   0.779  0.43643    
## stat_y1.l4 -0.0101087  0.0052842  -1.913  0.05663 .  
## stat_z1.l5  0.0134695  0.0665529   0.202  0.83974    
## stat_x1.l5 -0.0011361  0.0302069  -0.038  0.97002    
## stat_y1.l5 -0.0028914  0.0052850  -0.547  0.58470    
## stat_z1.l6  0.0039481  0.0664948   0.059  0.95269    
## stat_x1.l6 -0.0260481  0.0304551  -0.855  0.39302    
## stat_y1.l6  0.0033576  0.0052542   0.639  0.52326    
## stat_z1.l7 -0.1326844  0.0666843  -1.990  0.04746 *  
## stat_x1.l7  0.1110872  0.0301637   3.683  0.00027 ***
## stat_y1.l7 -0.0009127  0.0051652  -0.177  0.85985    
## stat_z1.l8  0.1557627  0.0568333   2.741  0.00647 ** 
## stat_x1.l8 -0.0441047  0.0285825  -1.543  0.12379    
## stat_y1.l8 -0.0103884  0.0045312  -2.293  0.02251 *  
## const       0.0003597  0.0025289   0.142  0.88697    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.04715 on 325 degrees of freedom
## Multiple R-Squared: 0.4727,  Adjusted R-squared: 0.4338 
## F-statistic: 12.14 on 24 and 325 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation stat_x1: 
## ======================================== 
## stat_x1 = stat_z1.l1 + stat_x1.l1 + stat_y1.l1 + stat_z1.l2 + stat_x1.l2 + stat_y1.l2 + stat_z1.l3 + stat_x1.l3 + stat_y1.l3 + stat_z1.l4 + stat_x1.l4 + stat_y1.l4 + stat_z1.l5 + stat_x1.l5 + stat_y1.l5 + stat_z1.l6 + stat_x1.l6 + stat_y1.l6 + stat_z1.l7 + stat_x1.l7 + stat_y1.l7 + stat_z1.l8 + stat_x1.l8 + stat_y1.l8 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## stat_z1.l1  0.0380076  0.1192246   0.319  0.75009    
## stat_x1.l1  0.5428587  0.0580590   9.350  < 2e-16 ***
## stat_y1.l1  0.0023643  0.0093319   0.253  0.80015    
## stat_z1.l2  0.4110311  0.1376729   2.986  0.00305 ** 
## stat_x1.l2 -0.1239843  0.0621466  -1.995  0.04687 *  
## stat_y1.l2 -0.0076661  0.0107074  -0.716  0.47453    
## stat_z1.l3 -0.2177210  0.1394022  -1.562  0.11930    
## stat_x1.l3  0.1559805  0.0622659   2.505  0.01273 *  
## stat_y1.l3  0.0113732  0.0109301   1.041  0.29886    
## stat_z1.l4  0.2245877  0.1394534   1.610  0.10826    
## stat_x1.l4  0.0111100  0.0625426   0.178  0.85912    
## stat_y1.l4 -0.0112230  0.0109834  -1.022  0.30763    
## stat_z1.l5  0.0327797  0.1383319   0.237  0.81283    
## stat_x1.l5 -0.0757495  0.0627859  -1.206  0.22851    
## stat_y1.l5  0.0106118  0.0109851   0.966  0.33475    
## stat_z1.l6 -0.1475144  0.1382112  -1.067  0.28662    
## stat_x1.l6 -0.0321001  0.0633017  -0.507  0.61243    
## stat_y1.l6  0.0067477  0.0109209   0.618  0.53709    
## stat_z1.l7  0.2933078  0.1386050   2.116  0.03509 *  
## stat_x1.l7  0.4440817  0.0626961   7.083 8.78e-12 ***
## stat_y1.l7  0.0009451  0.0107361   0.088  0.92991    
## stat_z1.l8 -0.1966949  0.1181295  -1.665  0.09686 .  
## stat_x1.l8 -0.3709662  0.0594095  -6.244 1.33e-09 ***
## stat_y1.l8  0.0069602  0.0094182   0.739  0.46043    
## const      -0.0002448  0.0052563  -0.047  0.96288    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.09799 on 325 degrees of freedom
## Multiple R-Squared: 0.5177,  Adjusted R-squared: 0.4821 
## F-statistic: 14.54 on 24 and 325 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation stat_y1: 
## ======================================== 
## stat_y1 = stat_z1.l1 + stat_x1.l1 + stat_y1.l1 + stat_z1.l2 + stat_x1.l2 + stat_y1.l2 + stat_z1.l3 + stat_x1.l3 + stat_y1.l3 + stat_z1.l4 + stat_x1.l4 + stat_y1.l4 + stat_z1.l5 + stat_x1.l5 + stat_y1.l5 + stat_z1.l6 + stat_x1.l6 + stat_y1.l6 + stat_z1.l7 + stat_x1.l7 + stat_y1.l7 + stat_z1.l8 + stat_x1.l8 + stat_y1.l8 + const 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## stat_z1.l1 -0.04632    0.78113  -0.059 0.952755    
## stat_x1.l1  0.36761    0.38039   0.966 0.334565    
## stat_y1.l1  0.62507    0.06114  10.223  < 2e-16 ***
## stat_z1.l2  2.21150    0.90200   2.452 0.014741 *  
## stat_x1.l2 -0.15593    0.40717  -0.383 0.701992    
## stat_y1.l2 -0.26992    0.07015  -3.848 0.000144 ***
## stat_z1.l3 -2.12668    0.91333  -2.328 0.020499 *  
## stat_x1.l3  0.83255    0.40795   2.041 0.042078 *  
## stat_y1.l3  0.16771    0.07161   2.342 0.019791 *  
## stat_z1.l4  0.62343    0.91367   0.682 0.495508    
## stat_x1.l4 -0.36084    0.40976  -0.881 0.379181    
## stat_y1.l4 -0.04551    0.07196  -0.632 0.527529    
## stat_z1.l5  0.33351    0.90632   0.368 0.713122    
## stat_x1.l5  0.11715    0.41136   0.285 0.775980    
## stat_y1.l5 -0.03512    0.07197  -0.488 0.625900    
## stat_z1.l6  0.54496    0.90553   0.602 0.547720    
## stat_x1.l6 -0.47645    0.41474  -1.149 0.251487    
## stat_y1.l6  0.07253    0.07155   1.014 0.311492    
## stat_z1.l7 -0.22129    0.90811  -0.244 0.807633    
## stat_x1.l7  1.57219    0.41077   3.827 0.000155 ***
## stat_y1.l7  0.28227    0.07034   4.013 7.45e-05 ***
## stat_z1.l8 -0.65807    0.77396  -0.850 0.395800    
## stat_x1.l8 -0.83627    0.38924  -2.148 0.032415 *  
## stat_y1.l8 -0.19753    0.06171  -3.201 0.001504 ** 
## const      -0.01403    0.03444  -0.407 0.684052    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.642 on 325 degrees of freedom
## Multiple R-Squared: 0.5114,  Adjusted R-squared: 0.4753 
## F-statistic: 14.17 on 24 and 325 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           stat_z1   stat_x1  stat_y1
## stat_z1 0.0022227 0.0009752 0.008576
## stat_x1 0.0009752 0.0096028 0.026390
## stat_y1 0.0085760 0.0263899 0.412208
## 
## Correlation matrix of residuals:
##         stat_z1 stat_x1 stat_y1
## stat_z1  1.0000  0.2111  0.2833
## stat_x1  0.2111  1.0000  0.4194
## stat_y1  0.2833  0.4194  1.0000
var = restrict(var_xx)

summary(var)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: stat_z1, stat_x1, stat_y1 
## Deterministic variables: const 
## Sample size: 350 
## Log Likelihood: 618.105 
## Roots of the characteristic polynomial:
## 0.9468 0.9468 0.9033 0.9033 0.8963 0.8963 0.8924 0.8924 0.8775 0.8775 0.8447 0.8447 0.8169 0.8169 0.7986 0.7986 0.7916 0.7916 0.7853 0.7853 0.7759 0.6891 0.6891 0.4986
## Call:
## VAR(y = XXX, lag.max = 8)
## 
## 
## Estimation results for equation stat_z1: 
## ======================================== 
## stat_z1 = stat_z1.l1 + stat_y1.l2 + stat_y1.l3 + stat_y1.l4 + stat_z1.l7 + stat_x1.l7 + stat_z1.l8 + stat_y1.l8 
## 
##             Estimate Std. Error t value Pr(>|t|)    
## stat_z1.l1  0.605130   0.040448  14.961  < 2e-16 ***
## stat_y1.l2 -0.008347   0.003476  -2.401  0.01686 *  
## stat_y1.l3  0.010248   0.003914   2.618  0.00923 ** 
## stat_y1.l4 -0.010841   0.003402  -3.187  0.00157 ** 
## stat_z1.l7 -0.116563   0.053594  -2.175  0.03032 *  
## stat_x1.l7  0.093560   0.021269   4.399 1.45e-05 ***
## stat_z1.l8  0.145807   0.053462   2.727  0.00671 ** 
## stat_y1.l8 -0.014068   0.003216  -4.374 1.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.04661 on 342 degrees of freedom
## Multiple R-Squared: 0.4582,  Adjusted R-squared: 0.4456 
## F-statistic: 36.16 on 8 and 342 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation stat_x1: 
## ======================================== 
## stat_x1 = stat_x1.l1 + stat_z1.l2 + stat_x1.l2 + stat_x1.l3 + stat_z1.l7 + stat_x1.l7 + stat_z1.l8 + stat_x1.l8 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## stat_x1.l1  0.55284    0.04785  11.554  < 2e-16 ***
## stat_z1.l2  0.33424    0.08844   3.779 0.000185 ***
## stat_x1.l2 -0.14326    0.05051  -2.836 0.004836 ** 
## stat_x1.l3  0.16129    0.04629   3.484 0.000558 ***
## stat_z1.l7  0.25198    0.11130   2.264 0.024207 *  
## stat_x1.l7  0.43684    0.04823   9.058  < 2e-16 ***
## stat_z1.l8 -0.22368    0.11107  -2.014 0.044815 *  
## stat_x1.l8 -0.33751    0.05031  -6.709 8.16e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.09707 on 342 degrees of freedom
## Multiple R-Squared: 0.5022,  Adjusted R-squared: 0.4905 
## F-statistic: 43.13 on 8 and 342 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation stat_y1: 
## ======================================== 
## stat_y1 = stat_y1.l1 + stat_z1.l2 + stat_y1.l2 + stat_y1.l3 + stat_x1.l7 + stat_y1.l7 + stat_x1.l8 + stat_y1.l8 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## stat_y1.l1  0.62692    0.05035  12.452  < 2e-16 ***
## stat_z1.l2  1.43554    0.59423   2.416 0.016225 *  
## stat_y1.l2 -0.24272    0.05525  -4.393 1.49e-05 ***
## stat_y1.l3  0.15000    0.04731   3.171 0.001658 ** 
## stat_x1.l7  1.34904    0.35491   3.801 0.000171 ***
## stat_y1.l7  0.32573    0.05543   5.876 9.95e-09 ***
## stat_x1.l8 -0.85976    0.35904  -2.395 0.017175 *  
## stat_y1.l8 -0.21285    0.05633  -3.779 0.000186 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.6398 on 342 degrees of freedom
## Multiple R-Squared: 0.4895,  Adjusted R-squared: 0.4775 
## F-statistic: 40.99 on 8 and 342 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##          stat_z1  stat_x1  stat_y1
## stat_z1 0.002286 0.001003 0.009026
## stat_x1 0.001003 0.009915 0.027122
## stat_y1 0.009026 0.027122 0.430602
## 
## Correlation matrix of residuals:
##         stat_z1 stat_x1 stat_y1
## stat_z1  1.0000  0.2106  0.2877
## stat_x1  0.2106  1.0000  0.4151
## stat_y1  0.2877  0.4151  1.0000

Test af residualerne

roots(var)
##  [1] 0.9468395 0.9468395 0.9033093 0.9033093 0.8963237 0.8963237 0.8924446
##  [8] 0.8924446 0.8774510 0.8774510 0.8446970 0.8446970 0.8169476 0.8169476
## [15] 0.7986099 0.7986099 0.7915700 0.7915700 0.7852532 0.7852532 0.7759217
## [22] 0.6891159 0.6891159 0.4985940
serial.test(var, type = "BG")
## 
##  Breusch-Godfrey LM test
## 
## data:  Residuals of VAR object var
## Chi-squared = 43.21, df = 45, p-value = 0.548
arch.test(var)
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object var
## Chi-squared = 203.1, df = 180, p-value = 0.1143
normality.test(var)
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object var
## Chi-squared = 177.88, df = 6, p-value < 2.2e-16
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object var
## Chi-squared = 43.454, df = 3, p-value = 1.971e-09
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object var
## Chi-squared = 134.42, df = 3, p-value < 2.2e-16

Impulse Response Funktioner

irf = irf(var, ortho = TRUE)
plot(irf)

Opstilling af SVAR

amat = diag(3)
diag(amat) = 1
amat[2,1] = NA
amat[3,1] = NA
amat[3,2] = NA

amat
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]   NA    1    0
## [3,]   NA   NA    1
svar = SVAR(var, Amat= amat, estmethod = "direct")

summary(svar)
## 
## SVAR Estimation Results:
## ======================== 
## 
## Call:
## SVAR(x = var, estmethod = "direct", Amat = amat)
## 
## Type: A-model 
## Sample size: 350 
## Log Likelihood: -1489.885 
## Method: direct 
## Number of iterations: 196 
## Convergence code: 0 
## 
## LR overidentification test:
## 
##  LR overidentification
## 
## data:  XXX
## Chi^2 = 4192, df = 3, p-value <2e-16
## 
## 
## Estimated A matrix:
##         stat_z1 stat_x1 stat_y1
## stat_z1  1.0000   0.000       0
## stat_x1 -0.4376   1.000       0
## stat_y1 -2.8767  -2.444       1
## 
## Estimated B matrix:
##         stat_z1 stat_x1 stat_y1
## stat_z1       1       0       0
## stat_x1       0       1       0
## stat_y1       0       0       1
## 
## Covariance matrix of reduced form residuals (*100):
##         stat_z1 stat_x1 stat_y1
## stat_z1  100.00   43.76   394.6
## stat_x1   43.76  119.15   417.1
## stat_y1  394.62  417.09  2254.6